!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! 
!! Turbulence closures: base layout
!!
!! \n
!!
!! See \c mod_turb_flux for the documentation.
!<----------------------------------------------------------------------
module mod_turb_flux_base

!-----------------------------------------------------------------------

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_messages, only: &
   mod_messages_initialized, &
   error, warning, info

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_atm_refstate, only: &
   mod_atm_refstate_initialized, &
   t_atm_refstate

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_turb_flux_base_constructor, &
   mod_turb_flux_base_destructor,  &
   mod_turb_flux_base_initialized, &
   c_turbmod, c_turbmod_progs, c_turbmod_diags, &
   t_turb_diags

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> Turbulence model
 type, abstract :: c_turbmod
  integer :: d !< dimension, stored here for convenience
  integer :: ntrcs !< number of tracers (stored for convenience)
  logical :: initialized = .false.
 contains 
  !> General allocations
  procedure(i_init),                deferred, pass(tm) :: init
  !> Clean-up all the objects allocated by \c init, except for the
  !! prognotic variables, since there can be many of them.
  procedure(i_clean),               deferred, pass(tm) :: clean
  !> Diagnose the \c ngrad variables used in \c grad
  !!
  !! \note All the input/output is intended at physical points (i.e.
  !! <em>not</em> in modal representation): this is very important
  !! because this function can be called for both sides and elements.
  !!
  !! \note Velocity, temperature and the remaining variables which we
  !! need the gradient of are not polynomial functions, since they are
  !! not the primal unknown of the formulation. In practice, we
  !! compute the gradient of the \f$L^2\f$ projections of such
  !! variables.
  procedure(i_compute_grad_diags),  deferred, pass(tm) ::              &
                                                      compute_grad_diags
  !> Update the turbulence coefficients and other diagnostic quantities
  procedure(i_compute_coeff_diags), deferred, pass(tm) ::              &
                                                     compute_coeff_diags
  !> Compute the turbulent flux
  !!
  !! This is the main function for the turbulent module, and can be
  !! also used to compute diagnostics.
  !!
  !! \note This function is always evaluated at a pair
  !! (element,coordinates), so that it can be used both inside an
  !! element and on the element boundary. The evaluation on the
  !! element boundary, in any case, is done from one of the element
  !! connected to the boundary, meaning that there is no treatment of
  !! the numerical fluxes in this subroutine. The numerical fluxes are
  !! computed in \c mod_dgcomp_rhs combining the values provided
  !! separately by this subroutine for the two elements connected to
  !! the side.
  procedure(i_flux),                deferred, pass(tm) :: flux
 end type c_turbmod

 !> Turbulence model, prognostic variables
 !!
 !! Since these variables solve a prognostic equation, all the methods
 !! defined for \c c_stv, required by the time integrators, must be
 !! defined.
 type, extends(c_stv), abstract :: c_turbmod_progs
 end type c_turbmod_progs

 !> Turbulence model, diagnostic variables
 !!
 !! Contrary to \c c_turbmod_progs, for these variables we don't need
 !! the vector space operations. Also, this type is not abstract,
 !! since the basic layout could suffice for the simples models.
 type :: c_turbmod_diags
  integer :: ngrad !< number of scalar gradients to be computed
  !> Gradient to be diagnosed
  !!
  !! Indexes are: spatial dimension, variables (\c ngrad), basis
  !! functions (<tt>base\%pk</tt>), and elements.
  !!
  !! This variable is defined here, since it is used by the turbulence
  !! model, but is computed in mod_dgcomp_rhs, because the whole
  !! finite element concept is required.
  real(wp), allocatable :: grad(:,:,:,:)
 end type c_turbmod_diags

 !> Turbulence diagnostics
 !!
 !! Collect the turbulence field diagnostics. This module in not
 !! included in \c c_turbmod to allow specification of different
 !! intents for the turbulence model and the diagnostics.
 !!
 !! \note The diagnostics are not, in general, polynomial functions,
 !! so that the representation in terms of the finite element basis is
 !! not exact. In particular, it is possible that the finite element
 !! representation displays negative values also for diagnostics that
 !! should be nonnegative by definition, and such negative values
 !! could appear also at the quadrature points. However, this is not a
 !! problem for the code, which never uses the finite element
 !! representation of the coefficients, but only the coefficient
 !! values at the quadrature nodes.
 type :: t_turb_diags
  integer :: ndiags !< number of diagnostics
  character(len=100), allocatable :: diag_names(:)
  !> Modal values for the defined diagnostics
  real(wp), allocatable :: diags(:,:,:)
 end type t_turb_diags


 abstract interface
  pure subroutine i_init(tm,progs,diags,td,grid,base,ntrcs)
   import :: t_grid, t_base, t_turb_diags, c_turbmod, &
     c_turbmod_progs, c_turbmod_diags
   implicit none
   type(t_grid),       intent(in) :: grid
   type(t_base),       intent(in) :: base
   integer,            intent(in) :: ntrcs
   class(c_turbmod),   intent(inout) :: tm
   class(c_turbmod_progs), allocatable, intent(out) :: progs
   class(c_turbmod_diags), allocatable, intent(out) :: diags
   type(t_turb_diags), intent(out) :: td
    ! A typical implementation, for a model which needs the gradients
    ! of temperature and velocity (plus tracers), would be
    !
    ! tm%d           = grid%d
    ! tm%ntrcs       = ntrcs
    ! tm%initialized = .true.
    !
    ! allocate(selected_tm_diags::diags)
    ! diags%ngrad    = 1 + tm%d + tm%ntrcs
    ! allocate(diags%grad(tm%d,diags%ngrad,base%pk,grid%ne))
    !
    ! if(this_model_requires_progs_variables) &
    !   allocate(selected_tm_progs::progs)
  end subroutine i_init
 end interface

 abstract interface
  pure subroutine i_clean(tm,diags,td)
   import :: t_turb_diags, c_turbmod, c_turbmod_progs, c_turbmod_diags
   implicit none
   class(c_turbmod),   intent(inout) :: tm
   class(c_turbmod_diags), allocatable, intent(inout) :: diags
   type(t_turb_diags), intent(out) :: td
  end subroutine i_clean
 end interface

 abstract interface
  pure subroutine i_compute_grad_diags(tm,gd , consv,atm_ref)
   import :: wp, t_atm_refstate, c_turbmod
   implicit none
   class(c_turbmod),      intent(in)  :: tm
   class(t_atm_refstate), intent(in)  :: atm_ref(:)
   !> conservative variables (deviations)
   real(wp),              intent(in)  :: consv(:,:)
   real(wp),              intent(out) :: gd(:,:) !< grad variables
   ! \todo in general, one should probably add an input argument with
   ! the additional prgnostic variables; however, be careful how the
   ! input values are selected: element/side dofs/values?
  end subroutine i_compute_grad_diags
 end interface

 !> Diagnostic variables, excluding the gradients which are computed
 !! in \c mod_dgcomp_rhs with calls to compute_grad_diags. Notice that
 !! \c cd has input inout, since the \c grad field typically is
 !! already defined when this function is called, and must be
 !! preserved.
 abstract interface
  pure subroutine i_compute_coeff_diags(tm,cd , grid,base, &
                                        uuu,progs,atm_ref)
   import :: c_turbmod,c_turbmod_diags, t_grid,t_base, &
     wp,c_turbmod_progs,t_atm_refstate
   implicit none
   class(c_turbmod),       intent(in)    :: tm
   type(t_grid),           intent(in)    :: grid
   type(t_base),           intent(in)    :: base   
   real(wp),               intent(in)    :: uuu(:,:,:)
   class(t_atm_refstate),  intent(in)    :: atm_ref(:,:)
   !> This argument is allocatable for the following reason: passing
   !! unallocated allocatable is only allowed if the dummy argument is
   !! also allocatable, hence, without the allocatable attribute, one
   !! would have to define an allocated actual argument, which implies
   !! extending \c c_turbmod_progs. With the present implementation,
   !! when there are no prognostic variables, one can simply pass an
   !! unallocated object with the abstract type. Basically, the
   !! allocatable attribute has a role comparable to "optional".
   !! Notice that this is one of the very few cases when there is a
   !! good reason for having an allocatable, intent(in) dummy
   !! argument.
   class(c_turbmod_progs), allocatable, intent(in) :: progs
   class(c_turbmod_diags), intent(inout) :: cd
  end subroutine i_compute_coeff_diags
 end interface

 !> This interface includes some variables which are typically
 !! available during the computation and thus should not be
 !! recomputed: \c rho, \c p, \c uu, \c cc.
 !!
 !! \warning If the optional argument \c td is present, it is assumed
 !! that the function is called for all the quadrature points in one
 !! element (i.e. not for a side).
 abstract interface
  pure subroutine i_flux(fem, tm,ie,x,bp, consv,atm_ref, rho,p,uu,cc, &
                         progs,diags, td)
   import :: wp, c_turbmod, t_atm_refstate, t_turb_diags, &
     c_turbmod_progs, c_turbmod_diags
   implicit none
   class(c_turbmod),      intent(in) :: tm
   integer,               intent(in) :: ie
   real(wp),              intent(in) :: x(:,:)  !< coordinates
   real(wp),              intent(in) :: bp(:,:) !< basis functions
   !> conservative variables (deviations)
   real(wp),              intent(in) :: consv(:,:)
   class(t_atm_refstate), intent(in) :: atm_ref(:)
   !> Total values for density, pressure, velocty and tracers
   real(wp),              intent(in) :: rho(:), p(:), uu(:,:), cc(:,:)
   !> See the comments in \c i_compute_coeff_diags concerning the
   !! allocatable attribute.
   class(c_turbmod_progs), allocatable, intent(in) :: progs
   class(c_turbmod_diags), intent(in) :: diags
   real(wp),              intent(out) :: fem(:,:,:)
   !> Compute diagnostics
   type(t_turb_diags),    intent(inout), optional :: td
  end subroutine i_flux
 end interface
 
! Module variables

 ! public members
 logical, protected ::               &
   mod_turb_flux_base_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_turb_flux_base'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_turb_flux_base_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
  (mod_state_vars_initialized.eqv..false.) .or. &
        (mod_grid_initialized.eqv..false.) .or. &
        (mod_base_initialized.eqv..false.) .or. &
(mod_atm_refstate_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_turb_flux_base_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_turb_flux_base_initialized = .true.
 end subroutine mod_turb_flux_base_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_turb_flux_base_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_turb_flux_base_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_turb_flux_base_initialized = .false.
 end subroutine mod_turb_flux_base_destructor

!-----------------------------------------------------------------------

end module mod_turb_flux_base

